Institute of Physics, Jagiellonian University, 
Reymonta 4, 30-059 Krakow, Poland 



^ ■ ON THE PROPAGATION OF NON STATIONARY 

PRESSURE WAVES IN STELLAR INTERIORS 
(N 

"tj ■ Patryk Mach" 1 " 

o 

(N 

> 

j-^v ' An analysis of the propagation of non stationary waves in the adia- 

^f-\ , batic region of stellar interior is presented. An equation of motion with 

\Q • an effective potential is derived, similar to the Zerilli equation known in 

C^) ' the propagation of gravitational waves. The Huyghens principle is violated 

in this case and the energy diffusion outward null cones is expected. Nu- 
merical calculations demonstrate that the diffusion is weak for the case of 
C"| . standard Solar model; thus no significant effect corresponding to quasinor- 

O ^ mal modes can be expected. The likely reason for the absence of stronger 

features is the restriction of our analysis to adiabatic regions only, where 
i, • the breakdown of the Huyghens principle is insignificant. 

CO 

1. Motivation 

S^ . In the standard helioseismology (or asteroseismology, to be more gen- 

eral) one usually deals with stationary perturbations of the gas pressure 
(or density) in stellar interiors. This issue has been well investigated so 
far, both from the theoretical and observational side (see e.g. {2]). Here, 
by stationary approach we understand posing a hydrodynamical boundary 
problem leading to some characteristic frequencies that can be subsequently 
compared with the frequencies of observed stellar oscillations. There is no 
reason, however, why not to consider the propagation of non stationary 
waves in stellar interiors. It is known that waves propagating in an inho- 
mogeneous medium can produce some non stationary effects such as, for 
instance, appearance of the quasinormal modes (for an example taken from 
the theory of perturbations of the Schwarzschild space-time see jS])- These 
are especially important from the observational point of view as their fre- 
quencies and damping coefficients are independent of the wave profile but 



e-mail: mach@th.if.uj.edu.pl 

(1) 



depend on the characteristics of the medium. It is also known that in some 
cases the quasinormal modes can dominate (!]. 

In this paper we perform a simplified analysis of the problem which 
appears in the astrophysics of stellar interiors. The order of this work is 
as follows. In Section 2 we recall some basic formalism commonly used 
in the theory of stellar oscillations. Section 3 gives the description of the 
propagation of non stationary waves in the adiabatic region of the stellar 
interior together with the derivation of the exact form of the equation of 
motion. We transform this equation to the form of the Zerilli equation [7j, 
[S] which one encounters in the theory of gravitational waves propagating 
in the Schwarzschild background metric (a case known to give significant 
quasinormal modes). The Lagrangian formulation of the problem is given 
in Section 4 while in Section 5 we deal with the Noether's energy density 
and its diffusion through the characteristics. In Section 6 we present the 
effective potential occurring in the equation of motion, obtained for the 
case of standard solar model. Finally Section 7 shows the results of some 
numerical calculation of the propagation of non stationary pressure waves 
in the Sun. Some final remarks and conclusions are given in Section 8. 

2. The formalism 

In this section we shall remind some basic equations known from the 
Newtonian theory of stellar oscillations (or asteroseismology) . We will not 
give any precise derivation here as it can be easily found in other papers. 
In turn, we will try to focus our attention mainly on putting down all the 
assumptions leading to the mentioned equations and we will explain all the 
notation we use. 

We will consider the motion of gas in the star ruled by the Euler equation 

gd t v + £>vVv = -Vp + gg (1) 

together with the continuity equation 

d t Q + V(-v 6 ) = 0. (2) 

Here g denotes the density, p the pressure and v the velocity field of the gas. 
The term g in the equation (^Q) stands for the gravitational acceleration. In 
further considerations it will be, however, much more convenient to use the 
gravitational potential $ instead of g. We shall assume that a non minus 
convention g = V<3? holds. Then, the potential $ satisfies Poisson's equation 
of the form 

V 2 $ = -AttGq. (3) 



The above set of equations should be completed with one more, namely 
energy conservation equation (or the first law of thermodynamics) 

dq 1 f dp Fipdg^\ 



dt giX'i ~ 1) \ dt g dt 

Here q denotes specific heat (i.e. heat per unit mass) and we have used 
standard thermodynamic notation (see e.g. (3]) 
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The next step is to derive a set of linearized equations describing the 
evolution of small perturbations of the equilibrium structure of the star. 
Under assumptions of adiabaticity of the motion and the spherical symmetry 
of the undisturbed medium one may obtain: 

g d^ = -d rP ' + g d r ^-g'g , (5) 



(6) 



Q = =; V + Qvi p ~ A -, — , (7) 

Ti,oPo V r i,o dr dr ) 

^d r (r 2 d r &) + V 2 h & = -4TrGg'. (8) 

The notation used here requires some explanation. The primed quantities 
denote Eulerian perturbations and are functions of both position r and time 
t, whereas quantities with zero indices describe the undisturbed medium 
(i.e. an equilibrium structure of the star) and are only functions of the 
distance from the center of the star, due to the spherical symmetry we have 
assumed. Thus, for instance 

g(r,t) = g (r) + g'(r,t), 

p(r,t) = p (r)+p'(r,t),... (9) 

By £ we have denoted the radial part of the gas displacement, i.e. 

6T = £er + Sh, (10) 

where e r stands for the unit vector in the radial direction and £/j is a hor- 
izontal part of the displacement vector. Finally V^ denotes the horizontal 
part of the gradient operator. The reader interested in rigorous derivation 
of the above equations may consult 0. 



3. Non stationary perturbations, effective potential 

The equations introduced in the preceding section possess a class of 
stationary solutions, which is actually one of the main interests of the theory 
of stellar oscillations. We may, however, try to look for the non stationary 
solutions that could, in fact, have some physical meaning. The aim of this 
section is to derive some simplified equations in the form suitable for further 
numerical search for such solutions. 

We proceed with the separation of variables. All perturbations of our 
interest, such as £, p', g' may be expanded in the series of spherical harmon- 
ics. To simplify the notation we will drop the adequate spherical harmonics 
indices in the amplitudes. Due to the linearity of the obtained equations it 
is sufficient to write 

p'{r,d,(j),t) = p(r,t)Y lm (0, $),... 
Taking into the account that 

o 1(1+1) 

vlY = ~^y,, 

and substituting the above expressions to the equations ©, ©, an d © 
we get 

SodU = -drP + Qod r ® - gg , (11) 

" d t6 ~ ^d r (r 2 Qo dtt) = ^4^(p ~ So*), (12) 

\d r (r 2 d r $) - ^-±iU = ~^ G ~e ( 13 ) 
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and 



P + QoUt; 1 1 — • (14) 



Ti,oPo V r i,o dr dr 

The mentioned stationary solutions can now be obtained by setting 

f(r,t) = f(r)e- iut 

for each amplitude f(r,t) of a hydrodynamical variable /. By completing 
ordinary differential equations obtained this way with the suitable boundary 
conditions we can determine characteristic frequencies u of the oscillation 
modes. 

We will, however, try to proceed in a different way. Instead of posing an 
eigenvalue, boundary problem, we will try to formulate some Cauchy prob- 
lem with just one, second order, partial differential equation, describing the 



time evolution of an initial perturbation. Additionally, we will not consider 
any boundary conditions and thus we will treat a star as a formally infinitely 
distributed medium. Of course, we will not consider the propagation of the 
perturbations beyond the assumed radius of the star. 

We shall now formulate two important, simplifying assumptions: 

1. We will put $ = everywhere. This assumption was examined for the 
first time by Cowling. It means simply neglecting the changes in the 
gravitational field that arise due to small perturbations of the density 
in the star. 

2. We will assume that the medium is adiabatic, i.e. 

1 d lnpo d In qq 
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dr dr 



0. (15) 



With this assumption we eliminate the propagation of the internal 
gravity waves, focusing attention only on the pressure waves in the 
star. Now, of course, we will have to examine, whether the star re- 
gion in which the perturbation propagates is indeed adiabatic. The 
convective zone in the stellar atmosphere can serve as an example of 
such region. 

Under above assumptions the considered set of equations reduces to the 
following form 

" #5 " ^dr(r 2 Q dtt) = l -^P, (16) 

Qodti= -d r p-Qgo, (17) 

Q0 p. (18) 
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We may now put the equation (|T8|) into the equation (|T6|) to obtain 

-d t p-- s d r (r QoOtt) - o P- 



i,oPo r r 2 

The term <9|£ in the last equation can be subsequently eliminated using (|17|) . 
Taking into account that —goQo = dpo/dr and remembering the condition 
(|15|) we get after some calculations 

-d^p + d?,p + [ j— ) d r p J 
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We can now get rid of the term containing the derivative d r p by introducing 
a new dynamical variable P, defined with the equation 



p = - P. 

r 

Indeed, after some calculations we obtain an equation of the form 

- -^-d 2 p + d r p - v ^ p = °> ( 19 ) 

c z (r) 
where 

c 2 (r) = hm. (20 ) 

By V we have denoted here a function playing the role of an effective po- 
tential 

V ( r \ _ l dln Qo l ^ 2 lngo 1 ( dlng \ 2 1(1 + 1) 
r dr 2 dr 2 4 \ dr J r 2 

Equation ()19|) can be still transformed into even more convenient form. We 
will eliminate the term c~ 2 standing before the time derivative of P by 
introducing a new coordinate 



dr' 

c(r') 
o 

Thus an obvious relation 

dr* 1 



(22) 



dr c(r) 
holds and the equation (fTHjl may be written as 

-d 2 t p + d 2 r *p - ^d r *p - c 2 vp = o. 

dr* 
We will get rid of the term proportional to d r * P in a way similar to that we 
had used before. We define a function II with a relation 

P = y/cR 
to derive the final version of our equation of motion 

- d 2 n + d 2 * n - vn = o (23) 

in which a new effective potential 

,„ 1 fd\nc\ 2 ld 2 \nc . t . 

V = ^ V + 4{^) -2*^ (M) 

has been introduced. 



4. Lagrangian description, energy 

One may notice that the equation (|23|) can be obtained from the varia- 
tional principle 

by taking an action S of the form 

S = J Cdtdr* = ~\ I ((d t II) 2 - (d r *U) 2 - VU 2 ) dtdr* . (25) 

The equation (|23|) appears then as the Euler-Lagrange equation for the 
Lagrangian density C, i.e. 

8 » £ - 8 'Mn- a ^ = - (26) 

We can now make use of the first Noether theorem applied to the action 
(|25|) . what should allow us to define an energy for the LT amplitudes. We will 
present this issue in a little detail due to some subtle matters that appear 
here. 

We begin with considering an infinitesimal time translation of some do- 
main £1 

^:I 2 D(]3 (t, r) h-> (t' , r) = (t - e, r) £ M 2 , 

which will be assumed to be a symmetry which means that a variation of the 
action (|25|) caused by this transformation vanishes in the domain £1. If in 
addition we assume that the motion happens to be real (the Euler-Lagrange 
equation (|26|) is satisfied) then after some calculations we obtain 

5S ~'L (( £ " 5i 8 ' n ) *" " (sfn 8 ' n ) *) - "■ < 27 » 

The amplitude II is defined in a half plane r* ^0. As a domain f2, over 
which we proceed with integration we may now take a part of that half plane 
enclosed between two constant time lines, given by the equations t = t\ and 
t = t 2 . 

Let us notice, that the amplitudes LT have appeared in the separation 
of variables in the (3 + 1) dimensional problem and, therefore, have to 
satisfy some additional conditions. In particular p needs to be finite in its 
domain and thus also at r* = r = 0. It follows simply that an equality 
II(r* = 0, t) = 0, and so d t Ii(r* = 0, t) = must hold. Therefore, for the Q 
chosen above we have 

dt= / — — -d t U dt = 0. 



dn \dd r *UJ J \dd r *U /r , =Q 



Finally it follows that the quantity 

oo 


is conserved, i.e. constant in time. The expression 



£ = £ - ^^n = - [[d t ny + (d r *ny + vu z j (28) 

can thus be interpreted as an energy density. 

In consistency with the comment made earlier, we have assumed here 
that a medium in which the waves propagate is infinite and all perturbations 
vanish at least at infinity. 

5. Energy diffusion 

Let us consider now 7, being a part of the incoming characteristic of 
the equation (|23j) . with an origin in the point with the coordinates r* = n, 
t = t\. The characteristic 7 divides the domain Vt enclosed between two 
constant time lines t = t\ and t = ti into two subdomains: the inner one, 
fil, and the outer one JI2 (Fig. ^). Let us consider next a point with the 
coordinate r* = R lying on the 7. The expression 

00 
j(R,t) = (-d R + d t ) fsdr* (29) 

R 

may be interpreted as a rate of the energy change along 7. A straightforward 
calculation making use of equation (|28|) and of the motion equation (|23j) 
shows that 

j(R,t) = U(d t u-d r ^u) 2 + vu 2 ) 

Z \ / r*=R 

Calculations leading to the above result may be simplified even more by 

noticing that 

00 

R 

what in fact we had obtained earlier by writing formula (|27|) . 

The whole energy AE, that diffused from a domain S7i to the domain 
O2 during a time between t\ and £2 may be calculated as an integral of the 




Fig. 1. An illustration of the scattering process considered in this paper. 



expression l)29JI over 7 



AE 






R=-t+r 1 +t 1 dt 



(30) 



t2 



. (dtu - d r *ny + vn 2 \ dt. 

2 \ /r*=-t+n+ti 



h 



Let us now imagine that an incoming perturbation of compact support 
enclosed initially in the domain f^i moves along the characteristic 7. The 
quantity AE given by the expression ()3(Jj) would represent the whole energy 
scattered outwards (into the domain Q2) during the period between t\ and 
£2- This scattering is mathematically due to the non- vanishing potential 
V and, to be more precise, to a potential that differs from a single term 
proportional to the 1(1 + l)/r 2 which will always arise from the separation 
of variables in the wave problem of spherical symmetry. 

Experience tells us that robust energy diffusion signals the presence of 
quasinormal modes. Since quasinormal modes are rather difficult to be 
found numerically, we prefer to deal with examining the energy diffusion 
instead @J. 

In the next sections of this paper we will concern ourselves with ex- 
amining such energy scattering for a realistic case, namely standard solar 
model. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 



Fig. 2. The r* coordinate plotted versus radius r. The r values are expressed in 
the radius of the Sun units, while r* is just arbitrarily normalized. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Fig. 3. The values of N 2 /go, N being the Brunt-Vaisala, frequency, obtained for 
the used solar model and plotted versus r*. 



6. Effective potential for a standard solar model 

We will now apply the results of the preceding sections to a standard 
solar model. We have already stated, that a fundamental assumption of 
our simplified model is that of adiabaticity of the region in which the waves 
propagate, and that validity of this assumption need to be carefully exam- 
ined. The standard solar model has been chosen because of the existence 
of the convective zone in which the condition (|L3|) is satisfied up to a high 
degree. 
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0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 



Fig. 4. The effective potential Vq for the used solar model plotted as a function of 
r*. 



All our numerical calculations were made on the basis of a standard 
solar model computed by Bahcall, Pinsonneault & Basin jj. The obtained 
relation between the variable r* and the radius r is shown on the Fig. |^1 It 
should be noticed here, that the r* variable was normalized in such a way 
that it changes between the values and 1 for the used data range. Next, 
the function 



AT2 

90 



1 d In po din go 



IY 



dr 



dr 



is plotted versus r* in Fig. El We have adopted the notation N 2 /go here, as 
N corresponds to the well known Brunt-Vaisala frequency. This plot shows, 
that we may regard the area with approximately r* gj 0.58 as satisfying our 
assumption of adiabaticity. The effective potential, i.e. V for I = is, in 
turn, plotted in Fig.0J Unfortunately the only one interesting feature of this 
potential, that is a clear peak at r* ~ 0.5, remains outside the adiabaticity 
area. Therefore we can not consider any effects caused by the existence of 
this peak as being physically meaningful. 

One remark should be made here concerning the figures presented in 
this section. It is not possible to differentiate the data presented by Bahcall 
et ah. in any straightforward way to obtain any smooth enough functions 
and therefore some smoothing procedure appears to be necessary. A slight 
modification of the so called Savitzky-Golay filter (see e.g. [5]) was used 
here to obtain presented result. 
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7. Example of some numerical calculations 

Finally it was possible to examine, how an initially incoming wave pack- 
age evolves. Fig. El presents an example result of the numerical experiment 
explained in Section 5. 

For simplicity, only a spherically symmetric perturbation i.e. the case 
with / = is considered here. Propagation of all other modes may be 
examined in exactly the same way by taking the potential V for an arbitrary 
I number. 

Here, the characteristic 7 was chosen to originate at r* = 0.845 and 
the initial {i.e. for t = 0), purely incoming perturbation was taken to be a 
function of the shape defined by 

n(r ., t = 0)= J.W(^). if'eM, 

[ 0, otherwise, 

thus being just one, bell like part of the squared sine function, centered on 
a compact support [a, b]. This is, of course, a continuous and differentiable 
function. 

We have also examined the case with an initial perturbation in the form 
of standard C°° class function of a compact support. Such function may be 
constructed in a well known way with an use of the exponential function. 
It has also a bell like shape but it gives lower FWHM to support length 
ratio. It appears that the squared sine function is much more useful for 
our purpose as, basically, the amount of scattered energy increases with an 
increase of the FWHM of the initial perturbation. 

In the presented case, a and b were given the values a = 0.73 and 
b = 0.845 which correspond to the initial data support lying entirely in 
the Q\ domain (see Fig. ^). The plot presented on Fig. [5] shows the time 
derivative of the energy that has diffused from the domain 0,i to the outside 
domain VL2 divided by the whole initial perturbation energy. The variable r* 
and the units of t (approximately 2500 s per unit) are defined in such a way 
that the sound speed expressed in these coordinates equals unity. Thus, 
looking at Fig. |S] it is easy to see how far could the initial perturbation 
arrive for a given time. Clearly, the large peak in the scattered energy 
corresponds to the mentioned bump in the effective potential which, as it 
has been already stated, cannot be considered in a convincing way as giving 
any results of physical meaning. It is, however, a good example of effects 
that may arise in the propagation of the non stationary waves due to the 
inhomogeneity of the medium. 

In our considerations we are of course restricted to the area where the 
adiabaticity assumption is satisfied. In fact, even in this area some diffusion 
of energy does occur but on a negligible scale. Fig. shows the energy that 
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Fig. 5. The time derivative of the energy scattered backwards. Here energy is 
expressed in the units of total energy of the initial perturbation, which is actually 
conserved, as described in section 4. 




Fig. 6. The energy E scatt scattered through the characteristic expressed as a func- 
tion of time t for t < 0.2. Here E to t denotes total initial energy of the wave 
package. 



has diffused through the characteristic as a function of time. These are in 
fact the same data which we have already plotted on the Fig.[5]but this time 
restricted to the times lower than 0.2 what corresponds to the propagation 
in the adiabatic zone. 
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8. Final remarks 

It is already well known that non stationary waves can carry interesting 
information about an inhomogeneous medium in which they propagate 0, 
EH j |H| [I] • In this paper we have examined the propagation of such waves 
in the Solar convective zone by looking at the energy diffusion process. It 
appears that the energy scattering occurs with rather negligible efficiency 
and, consequently, we expect the quasinormal modes to be absent. This 
may, however, follow from the fact that we have restricted ourselves to the, 
perhaps not interesting, simplified case of adiabatic media. It is possible that 
the investigation of the full model which takes into account all important 
physical aspects would lead to some positive results. It is also possible that 
positive results can be obtained by repeating the calculations presented in 
this paper for the models of some other stars. 
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